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ABSTRACT 

Aims. Radiative feedback from massive stars is a key process to understand how HII regions may enhance or inhibit star formation 
in pillars and globules at the interface with molecular clouds. We aim to contribute to model the interactions between ionization and 
gas clouds to better understand the processes at work. We study in detail the impact of modulations on the cloud-HII region interface 
and density modulations inside the cloud. 

Methods. We run three-dimensional hydrodynamical simulations based on Euler equations coupled with gravity using the 
HERACLES code. We implement a method to solve ionization/recombination equations and we take into account typical heating 
and cooling processes at work in the interstellar medium and due to ionization/recombination physics. 

Results. UV radiation creates a dense shell compressed between an ionization front and a shock ahead. Interface modulations produce 
a curved shock that collapses on itself leading to stable growing pillar- like structures. The narrower the initial interface modulation, 
the longer the resulting pillar. We interpret pillars resulting from density modulations in terms of the ability of these density modula- 
tions to curve the shock ahead the ionization front. 

Conclusions. The shock curvature is a key process to understand the formation of structures at the edge of HII regions. Interface and 
density modulations at the edge of the cloud have a direct impact on the morphology of the dense shell during its formation. Deeper 
in the cloud, structures have less influence due to the high densities reached by the shell during its expansion. 

Key words. Stars: formation - HII regions - ISM: structure - Methods: numerical 



1. Introduction 

Radiative feedback from massive stars might be an important 
process to explain the star formation rates on galactic scales. 
Its role in complexe structures like giant molecular clouds 
is sti l l a matter of debate (see iDale et af] 120071: iPrice et aP 
|2QTqI: iDale & Bonnellll20TTh . When the UV radiation from the 



massive objects photoionize the surrounding gas, a "bubble" 
of h ot ionized gas exp ands around the star: the HII region 
(see iPurcell et alJ |2QQ9[ for example). While star formation is 
inhibited inside the bubble, the small-scale compression at the 
edge of the HII region, due to its expansion, seems to form 
, elongated structures (pillars) and glob ules in which the star 
. formation activity seems enhanced (see iDeharveng et ani201Ql 
for example). Different models investigate this process. 
First studies of HII regions (e.g. StromgrenI 119391: 
■ lElme green &Ladal Il977h show how UV radiation leads to 
' the formation of an ionization front and of a shock ahead. The 
gas is compressed between them, forming a dense shell which 
may lead to fragmentation, gravitational collaps e and star forma- 
tion, t his is the collect and collapse scenario (lElmegreen & Ladal 
Il977h . 

An other one was proposed by iBertoldil (Il989l) . the radiation- 
driven implosion scenario. He looks at the photoevaporation 
of spherical neutral clouds and finds that the ionization front 
drives a shock into the cloud leading to the compression of 
the initial structure into a compact globule. This scenario 
has been studied in detail with numerical simulations for 



years (see iLefloch^ B. &Lazarefi Il994l: IWilliams et all 1200 ll: 
iKesse l-Devnet & Burkert 2003). RecentlvT lBisbas et al.l (l2009l) 
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and iGritschneder et al.i (i2009») looked at the implosion of 
isothermal spherical clouds with smoo thed particle hvdrody - 
namics (SPH) codes. Using a grid code, Mackey & LimI (1201 Ol) 
produce elongated structures from dense spherical clumps using 
the shadowing effects of these structures. 

In the last decade, the importance of the initial turbulence in 
the cloud has been studied with three-dim ensional simulations, 
at different scales. iMellema et al.l (l2006h present simulations 
of the formation of the HII region with a good agreement with 
observations. The interpla y between ionization and magnetic 
fields has been stud ied by iKrumholz et al.l (l2007l) and later by 
[Arthur et"an (l201ll) , in the context of HII region formation, 
finding that magnetic fields ten d to suppr e ss the small-scale 
fragmentation. On larger scales iDale et af] (l2007l) look at the 
impact of ionization feedback on the collapse of molecular 
clouds findin g a sli ght enhancement of star formation with ion- 
ization while lpale'& Bonnell ( 2 011) finds almost no impact. On 
smaller scale jLora et al.l (l2009l) find that the angular momentum 
of the resulting clumps is preferentially oriented pe rpendicular 
to the incident radiation. IGritschneder et al.l (1201 Ol) show that 
pillars arise preferentially at high turbulence and that the 
line-of- sight velocity structure of these pillars differs from a 
radiation driven implosion scenario. 



These HII regions are seen in a lot of massive molecular 
clouds and are the object of a l arge number of obser vation- 
nal studies: in Rosett e nebula (ISchneider et al.l |201Q|), Ml 6 
(lAndersen et al.l l2004l) . 30 Doradus (IWalborn et al.l l2002h . 
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Carina nebula d^mith et alJ 120001), Elephant Trunk nebula 



('Reach eta 
( Jiang et al 



]2004 . Trifid nebula (iLefloch et alJ l2002h. M17 
]2002h or also the Horsehead nebula (iBowler et alJ 



I2OO9I) . Spitzer observation s provide also a wide r ange of HII 
regions studied in detail by lDeharveng et al] (l2010l) . Pillars and 
globules are often seen and present density clumps which may 
lead to star formation. The different scenarii described above 
can be constrained thanks to these observations. 

The present study focuses on "simple" situations in order to 
highlight the key mechanisms at work in the interaction between 
a HII region and a cloud. We present the numerical methods 
needed for this study and then two different set ups, cloud-HII 
region interface modulation and density modulation inside the 
cloud and finally a study of the different stage of evolution of 
the resulting pillars. A following paper will investigate these sit- 
uations in a more complete set up which will include turbulence 
inside the cloud. 



2. Numerical methods 

We consider a molecular cloud impacted by the UV radiation 
of a OB cluster to study how structures develop at the inter- 
face between the resulting HII region and the cloud. Subsection 
12. II describes the method used to simulate gas h ydrodynamic in 
the rn olecular clou d wi th the HERACLES code (iGonzalez et all 
l2QQ7h . Subsection 12.21 describes the numerical method used to 
take in account the UV radiation from the OB cluster and the 
resulting ionization/recombination reactions. Thermal processes 
from these reactions and the heating and cooling rates used in 
the molecular cloud are described in subsect. 12.31 



2.1. Hydrodynamic 

Our simulations are performed with the HERACLES codeQ. It is 
a grid-based code using a second order Godunov scheme to solve 
Euler equations. These equations are given in Eq.[T]in a presence 
of a gravitationnal potential O, constrained by the Poisson equa- 
tion: AO = AnGp. 



dpW 



dp 
ot 



-h V • [pV V -h = -pVO , 
ot 

^+V.[(£ + /7)V] = -p(V-V(D) + A(p,r) 
ot 



(1) 



We use an ideal gas equation of state, so that E = p/(y- 1) + 
0.5pV^ with y = 5/3. A(p,T) corresponds to the heating and 
cooling processes detailed in subsect. [23l leading to a poly tropic 
gas of index around 0.7 at equilibrium (for p = 500 /crrP). All 
the other variables have their usual meaning. 

2.2. Ionization 

The stationnary equation governing radiative tranfer in spherical 
geometry is given by: 



Pdly/dr - P )/r X dly/dp = -Kyly -\- T]y 



(2) 



p is equal to cos^ and is the polar angle, Iv(r,p, v) is the 
spectral intensity at a radius r in the direction defined by p at the 



frequency v, Ky(r, v) is the absorption coefficient and 7/y(r, y) is 
the source term. Radiation is coming from a point source such 
that Iy(r,p, y) = Fy(r, v) X Si(p), in which Si(p) is a Dirac dis- 
tribution whose peak is at yu = 1 corresponding to the direction 
^ = in spherical coordinates. Integrating Eq.[2]over p leads to 
Eq.[3]in which the absorption coefficient and the source term are 
given by Eq.m 



dFyldr -\-2/rX Fy = -KyFy -\- rjy 



T]y(r,y) = S^(y)6o(r) 
Ky(r,y) = crriHyi/y^ 



(3) 



(4) 



hyi is equal to 13.6eV, a is the ionizing cross section for pho- 
tons at frequency vi and uhq^x) the density of neutral hydrogen 
atoms, ^o(^) is a Dirac peak at r=0 and is the flux from the 
central OB cluster. The radiation from the cluster is approxi- 
mated by an unique source of radius at a distance ro of the 
computational domain, corresponds to a diluted planckian 
function at the star temperature: S ^{y) = 7i{rllr^)By{T^). By di- 
viding Eq.[3]with hy to have an equation on the number of pho- 
tons rather than radiative energies and averaging y between y\ 
and -h 00 in order to keep track only of the ionizing photons, we 
get: 



l/r^xd(r^Fy)/dr: 



-HHoCTyFy -\- So(r)S , 



(5) 



Fy is now the number of ionizing photons per unit of 
surface and time arriving in the radial direction and d-y is the 
average cross-section over the planckian source at and S * the 
rate of emission of ionizing photons by the stars. We can get 
plane-parallel equations in the limit r ^ 00. 

We now compute the equations for photo-chemistry. We 
set X to be the fraction of ionization X = nH^/nn in which 

= ^//+ + ^//o, it follows the dynamic of the gas as an ad- 
vected quantity. The variation of protons is the number of in- 
coming photons which are going to interact minus the number 
of protons which will be used for recombination. The number of 
interacting photons in a volume dvceii is given by the number of 
incoming photons threw a surface dsin multiplied by the prob- 
ability of interaction given by the cross-section, dP = o-yUHodr 
leading to: 



diXnn) 
dn^. , 



dUy 



■drih 



Pnn+ne-dt = pX n^jdt 

Fy{r)dsindtdPI dvceii = o-yFynuiX - X)ojdt 



(6) 



http://irfu.cea.fr/Projets/Site_heracles/index.hmtl 



The dilution term co = (dsindr/dvceii) is equal to r^j^/r^^n in spher- 
ical coordinates and to 1 in the plan-parallel limit. is the radial 
position of the entrance surface for incoming photons dsin and 
fceii is the radial position of the centre of the volume dvceii- P 
gives the rate of recombination and is equal to 2 x I0"^^r"^-^^ 
cnr' I s in which T is the tempera ture of therm odynamic equilib- 
rium between all the species (^see lBlacklfTosTh . 



2.3. Thermal processes 

Thermal processes are taken in account by adding the heating 
and cooling rate A(p, T) in the equation of energy conservation 
in Eq.m 

In the ionized phase, we consider two processes which have a 
major impact on the thermodynamic of the gas, the photoelectric 
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heating due to the massive stars UV flux and the cooling due 
to recombination of electrons onto protons. The energy ey(T^) 
given by ionizing photons is the integrated value of Iy/hy(hy - 
hvi) between vi and -\-oo on the incoming spectrum /y. In the 
following we will assume a planckian distribution By(T^) with 
= 35000K which gives ey(T^) = leV. Therefore the heating 
rate is given by Eq.[71 The cooling term is the loss of the thermal 
energy of electrons used by recombination given by Eq. [H 



*H = dny.^Jdt X = aj(l - X)nHFyCryey 



£ = dnH,,JdtxhT/(y-l) ■ 



■-PX^lhTliy- 



1) 



(7) 
(8) 



At the equilibrium between ionization and recombina- 
tion, the recombination rate is equal to the ionization rate 
driy^^^ = driH,,^. Therefore when thermodynamic equilibrium 
is achieved in the ionized phase the temperature is given by 
(y - l)ey/kt which is equal to 7736 Kelvin. The corresponding 
isothermal curve in the pressure-density plane is drawn in Fig.[T] 
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Fig. 1: Thermodynamic equilibrium in the pressure-density 
plane for the highly and weakly ionized phases 



In the weakly ionized phase, we simulate the radiative heat- 
ing and cooling of the interstellar mediu m (ISM) by adding 
five major processes following the work of lAudit & Hennebelld 
(I2005h and lWolfire et"aD (ll995ll2QQ3h : 

- Photoelectric heating by UV radiation with a flux equal to 
Go/ 1.7 in which Go is the Habing's flux 

- cooling by atomic fine- structure lines of Cll 

- cooling by atomic fine- structure lines of Ol 

- cooling by H (hya) 

- cooling by electron recombination onto positively charged 
grains 

The UV flux used in this phase is an ambiant low flux, 
additional to the UV flux Fy coming from the massive stars, 
which is used in our ionization process described in subsect. 
12.21 This heating and cooling function is only valid for the dense 
cold and weakly ionized phase. Therefore these processes are 
weighted by 1 - X and contribute only in the weakly ionized 
phase. In this phase, the ionization fraction used for the thermal 
processes is a function of the temperature and is given by 
IWolfire et al.l ([2003.) (typically around lO""^). The thermody- 
namic equilibrium in the pressure-density plane is given in fig. [T] 



The transition at the ionization front is very sharp as it can 
be seen on the tests in subsect. 12.41 Therefore the fraction of the 
gas in between the two phases is small and is not dynamically 
significant. 

2.4. 1D spherical test: HII region expansion in spherical 
geometry 
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Fig. 2: Expansion of an HII region in a homogenous medium 
(n//=500 cm"^, T^^25 K, S* = 10^^ Yiomzing/^)- Case 1 is a simula- 
tion in which the gas temperature is instantaneously put at 7736 
K when the ionization fraction is bigger than 0.5. Case 2 takes 
into account the time evolution of temperature given by Eq. [7] 
andlU The analytic solution is computed from Eq.[TOlwhen the 
thermal equilibrium is reached instantaneously. The Stromgren 
radius is given by Eq. [9j 

We first start by testing our numerical algorithms on a sim- 
ple ID spherical test which can be compared to analytical so- 
lution. Starting from a region with nn = 500 cm"^ and T ^ 25 
K, we switch on a photon flux coming from a typically 04 star 
with 5* = 10^^ Jionizingl s. At the beginning the typical time for 
ionization/recombination is much shorter than the hydrodynam- 
ical typical timescale. Therefore there is a first phase of devel- 
opment in which the hydrodynamic is frozen and the medium is 
ionized in a small sphere around the source. The radius of this 
sphere is controll ed by the Stromgren formula given in Eq.[9](see 
IStromgrenll 19391) . Ionization stops when all photons are used to 
compensate recombinations inside the sphere. In our example 
this radius is = 2.38 pc. 



I \ 
\47ipnH) 



1/3 



(9) 



If the thermal equilibrium between ionization and recombi- 
nation is reached instantaneously, it is possible to compute an 
analytic expre ssion of the time evolution of the ionization front 
(ISpitzedll978h : 



Ri = Rs(l - e-^'"') 



1/3 



(10) 



The development of the ionized sphere before hydrodynamic 
starts to matter is shown in Fig.O Two cases are treated numer- 
ically, first the thermal equilibrium is reached instantaneously 
and the gas temperature is switched at 7736 K when the ioniza- 
tion fraction is bigger than 0.5. It allows us to compare our nu- 
merical solution to the analytic formula given in Eq.[TOl the error 
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Fig. 3: Density (top) and temperature (bottom) profiles of the 
expansion of an HII region in a homogenous medium (n//=500 
cm"^, T^25 K, S* = 10^^ Jionizingl^)- Solid lines correspond to an 
adiabatic cold phase and dashed lines to a cooled cold phase with 
the cooling described in Subsect. 12.31 



is less than 0.5%. Second we resolve numerically the time evolu- 
tion of the temperature given by the balance between Eq.[7]and[8l 
The ionization front is slowed down by the explicit treatement of 
cooling due to recombination. The recombination rate increases 
when the temperature decreases leading to a lower penetration 
of ionizing photons. 

After this first step the energy deposited by the ionizing 
photons creates a spherical region with a high pressure. This re- 
gion will therefore start expanding and pushing the surrounding 
interstellar matter, creating a shock front ahead the ionization 
front. This process creates an expanding ionized cavity and a 
shell of compressed gas. Typical density and temperature pro- 
files are given in Fig.[3l Solid lines corresponds to the adiabatic 
case for the cold gas and the dashed lines to the cooling function 
described in Subsect. [23l The shock front is cooled down by the 
thermal processes and the shell is compressed up to two order of 
magnitude. This phenomenon is at the origin of the idea of the 
collect and collapse sce nario for triggered star form ation around 
HII regions proposed by lElmegreen & Ladal (1 19771) . 

We can derive the equations governing the density, the 
speed and the thickness of the shell in the approximat ion of 
an iso t hermal D-criti cal shock following Elmegreen & Ladal 
(1 19771) : ISpitzej (Il978l) . The density and the pressure when the 
HII region has a radius r can be derived from the equilibrium be- 



tween ionization and recombination in the sphere and are given 
by: 



nil 
Pii 



noirs/rf 
2nuhTu 



(11) 



a is equal to 1/3 and to (3S ^/ 47111^3)^^^ in spherical geome- 
try while, in the plan-parallel limit, a is equal to 1/2 and to 
Fy/ri^. The parameters of the shell can be computed using Eq. 
nUand assuming a shell temperature Ty/;^// and the corre spond- 
ing sound speed c shell (for details, see lElmegreen & Ladall 19771) . 
They are given by Eq.[T2lin which n shell is the maximum density 
in the shell, Vsheii the shock speed, and Isheii the typical width of 
the shell. 



= ypiilc]j^J{jdmH) 



^shell 

Vshell = (8/?///3po) 

hhell = 



1/2 



^rly^]helll^]hell 



(12) 



Assuming a temperature in the shell of 10 kelvin and a ra- 
dius of 20 parsecs for the HII region we get a shell density of 
3.2xl0'^cm"^, a speed of 3.2 km/s and a width of 0.5 parsec. 
These values are in good agreement with the shell parameters 
we obtain in the ID simulation in Fig.[3l the maximum density 
when the Hll-region radius is 20 parsecs is 2.9xl0'^cm"^, the 
shell speed 4.4 km/s and the shell width 0.27 parsec. We can see 
that our code can follow with very good accuracy the ionization 
front and the subsequent evolution of the dense shell. We will 
therefore now turn our attention to more interesting 3D cases. 



3. Forming pillars 

The interstellar medium has a very complex structure with im- 
portant inhomogeneities and large density fluctuations. The pas- 
sage from ID to 3D brings a lot of new degrees of freedom. 
It allows us to study how the shell resulting from the collect 
and collapse scenario may be perturbed. The possibility to have 
localized density gradients has been ex tensively studied with 
the radiation driven implosion scenar io _(|Mackev&Limll2010l: 
iBisbas et"an 120091: iGritschneder et al.ll2009l) and more recently 

r et all " 



within a turbulent media (IGritschneder et al.i2010l) . In this work, 
we first want to focus on more simple and schematic situations 
in which the physical processes at work can be identified and 
studied more easily. We will therefore study two idealized cases. 
First we will look at the interaction of the ionization front with a 
medium of constant density having a modulated interface. Then 
we will consider a flat interface but the presence of overdense 
clumps in the medium. 

3.1. From Interface modulations to pillars 

We first study the interaction between the ionization front and 
an interface modulated by an axisymmetric sinus mode with 
constant heigh (amplitude of 0.5 parsec) and four diff'erent base 
widths (see Fig. |4] and Fig. O. The size of the computational 
domain is adapted to the typical observed length scale of pillars, 
4 X 2 X 2 pc^, at a resolution of 400 x 200 x 200 corresponding 
to spatial resolution of 1 x 10"^ pc. The boundary conditions 
are periodic in the directions perpendicular to the ionization 
propagation direction, imposed where the flux is incoming and 
free flo w at the opposite. Ionization processes introduced in 
subsect. 12.21 are taken in the plan-parallel limit. The snapshots 
around 500 ky show that the narrower the initial structures, the 
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Fig. 5: Density cuts of four simulations of an ionization front propagation with a modulated interface. The UV flux at the top of the 
box is 10^ photons per second per square centimetre, the box is 4x2x2 parsecs (half of the cut is shown). The dense gas is initially 
at 500 H/cm^ at a temperature around 25 K (thermal equilibrium from Fig.[T]). The (base width)/height ratios of the modulations are 
respectively 2 for case (a), 1.4 for case (b), 0.9 for case (c) and 0.5 for case (d). 
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Fig. 4: 3d snapshot of density at 500 ky of an ionization front 
propagation with a modulated interface. The UV flux at the top 
of the box is 10^ photons per second per square centimetre, the 
box is 4x2x2 parsecs. The dense gas is initially at 500 H/cm^ 
at a temperature around 25 K (thermal equilibrium from Fig. 
[B. The (base width)/height ratios (w/h ratios hereafter) of the 
modulation is w/h=0.9. 



longer the resulting pillars. With a (base width)/heigh ratio (w/h 
ratio hereafter) of 0.5 we get a structure whose size was almost 
multiplied by 3.5 in 500 ky. Besides the initial structures have 
less and less mass with decreasing widths and still manage to 
form longer pillars as can be seen on figure ©. High densities 
or high initial mass are not needed to form structures which are 
going to resist the ionization. The important factor here is how 
matter is distributed in space and interact with the propagating 
shock. Very little and low-density material well distributed in 
space can result in a long pillar-like structure. We will call these 
structures pillars hereafter, however a detailed comparison with 
observations is needed to see whether these structures emerging 
from idealized scenarii are a good approximation of the reality 
or not. 

To study the properties of the pillars, we monitor their sizes 
and masses in Fig. [6l The size gain identified previously is a 
continuous process in time at an almost constant speed which 
depends very strongly on the width/heigh ratio. Narrow initial 
modulation have a very fast growth while larger one grow 
slowly, if at all. The mass increase relative to the initial mass of 
the structure is almost the same in all simulations and does not 
seems to depend on the initial w/h ratio. It reaches a factor of 
5-6 in 500 ky. 

Vertical profiles of the mass in the pillars (see Fig. (T]) show 
that the mass gain identified in Fig. [6] is accumulated at the 
base. The mass of the head of the pillar slightly increases and 
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Fig. 6: Monitoring of the size (top) and the mass (bottom) of 
pillars identified in Fig. [5] The size is calculated by the diff'erence 
of the average position of the ionisation front and the its position 
at the centre of the y-z plane (head of the initial structure). The 
mass is calculated in a box between the two vertical positions 
defined to calculated the size and with a width defined by the 
width of the initial structure and is normalized by the mass in 
the initial modulation. 



then remains stable during the simulation while the central part 
connecting the head and the base is stretched and its mass de- 
creases. The profiles of the vertical speed show that the bases of 
the pillars have a roughly constant speed and the size variations 
are due to vertical speed diff'erence of the heads. The diff'erences 
between the sizes of the pillars can be directly deduced from the 
velocity diff'erences of the heads: a velocity diff'erence of 0.5 
km/s during 500 ky gives a spatial diff'erence of 0.25 pc which 
is approximatively the diff'erences observed between the simula- 
tions with a width/heigh ratio of 2 and the one with a ratio of 1 .4. 

A closer look at the simulations (see Fig. [S]) shows that 
the motions perpendicular to the propagation direction are 
shaping the structures. These motions are triggered by the 
initial curvature of the structure, two cases are distinguished, an 
initially convex zone and a initially concave zone. A shock front 
on a convex zone will trigger its collapse, it is the phenomenon 
which takes place in the radiation driven implosion scenario 
and here at the head of the initial structure. The curvature of the 
initial structure leads to lateral convergent shocks that collide 
on themselves. However the ionization of a concave zone is 
quite diff'erent, it will dig a hole in the medium, and the gas 
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Fig. 8: Zoom on the ionization of convex and concave zone in the w/h = 0.9 simulation between 22 and 132 ky. Arrows represent 
the velocity field projected in the x-y plane. 



is pushed away lateraly. The velocity fields in Fig. [8] illustrate 
these phenomena, on a initially convex zone the velocities point 
inwards and lead to the collapse of the central structure whereas 
on the initially concave zone at the base they point outwards and 
dig a hole in the media. These phenomena are at the origin of 
the size and mass increase. Indeed, the collapse of the structure 
takes more time if the structure is wide, therefore it will be 
accelerated down longer, explaining the diff'erence of vertical 
speed of the heads of the pillars seen in Fig. [T] Besides, the 
holes on each side of the pillar (see Fig. [8]) gather the gas at the 
base of the pillar, explaining the mass increase in Fig. [3 



Besides the mass increase, there is an important density gain. 
The mass histogram in the pressure-density plane in Fig.[9]shows 
that the compression due to the heating from UV radiation in- 
creases the pressure of one order of magnitude from the ini- 
tial state which is at a density of 500 cm^ and at a pressure of 
10"^^ erg/cm^. The gas is then distributed in two phases at equi- 
librium, hot-ionized on the green-dashed straight line and cold- 
dense-unionized on the blue dashed curve. The gas in transition 
between them is seen with the blue contour. AID simulation of 
the collect and collapse process can increase the density up to 
5x10"^ cm^, this limit is drawn with the red dashed line in Fig.[9j 
The mass at the base of the pillar, due to the holes from the initial 
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Fig. 7: Mass profiles (top) and vertical velocity profiles (bottom) 
at three different times of the pillar in the w/h = 0.9 simulation. 
These profiles are made within the boxes defined in Fig. [6] 



concave parts and from the collapse of the lateral shocks, leads to 
a density enhancement of 1-2 orders of magnitude. Although the 
structures in our simulations are still slightly under their Jeans 
lengths, this process has a better chance to trigger star forma- 
tion at the base of the pillars than the simple collect and collapse 
scenario. 



3.2. From density modulations to pillars 

Clouds have irregular shapes but they have also unhomo- 
geneities. Therefore, after this study of interface modulations, 
we consider density modulations that we call clumps hereafter 
(with no reference to an observational definition). We define 
them as regions of constant overdensity. Figure [TOl shows the 
diff'erent cases we studied. The three first cases are clumps in a 
homogeneous cold cloud, with a density contrast of respectively 
2, 3.5 and 5 and the last one corresponds to the radiation driven 
implosion scenario, the clump is "isolated" in a hot low-density 
medium so that the shock forms directly on the structure, the 
density in the clump is the same as the one with a density con- 
trast of 5. In all cases the clumps are at pressure equilibrium with 
the surrounding matter. 

In the isolated case, the shock forms at the surface of the 
clump and the shell is therefore initially curved by the shape of 
the clump with a w/h ratio of 2. This roughly corresponds to the 
wider case studied in sect. 13. 1[ Since the shock is curved, it will 
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Fig. 9: 2D plot of the mass fraction at a given density and pres- 
sure for the w/h = 0.9 simulation at the end of the simulation 
(554 ky). Blue to black contours are increasing mass fraction 
contours (blue: lE-6, green: lE-5, yellow: 5E-5, orange: lE-4, 
red: 5E-4, black: 2E-3). The dashed-red line is the maximum 
density achieved in a plan-parallel ID simulation : 4.8x10^ cm^. 
The total mass at a density above this limit is around 10 solar 
masses 



collapse on itself an d form a pillar-l i ke structure, this case was 
studied in detail by iMackey & LimI (l2010l) : iGritschneder et al.l 
( 2009). However, with a widht/heigh ratio of 2 the structure is 
quite small and look quite similar to the w/h=2 case of the pre- 
vious section. However, since it is isolated, the accumulation of 
mass process at the base identified in the previous section can 
not take place. 

In the three other cases, the shell is flat when it is formed, 
therefore the outcomes of the simulations are not so clear. When 
the clump is not dense enough (first case n}lciump/^^cioud='^), it 
will not curve the shell enough to have lateral colliding shocks. 
Therefore in this case, there is no head in the final structure 
which is very small (i.e comparable to the initial clump size). 
In the two other simulations (nllciump/^^cioud='^-^ and 5), the 
clumps curve the shell enough to trigger the collision of vertical 
shocks thereby forming an elongated pillar. This is very close 
to what we observed previously on interface modulations. 
Besides, we can identify on the final snapshots a base for the 
pillar- structures formed by lateral holes in the cloud and the 
associated accretion mass process discussed above. 

The importance of the curvature eff'ect can be emphasized 
by comparing the isolated case and the n}lciump/T^^cioud=^ case. 
In both cases the density of the clump is the same, however 
when the shell is formed flat on a homogeneous medium, 
the dense clump will resist the shock to form an elongated 
modulation on the shock surface. In the isolated case, the 
shock is formed instantaneously curved on the structure with a 
w/h ratio of 2. Furthermore, contrary to the isolated case, the 
structures will grow in mass and size because of the connection 
with the cloud. These efl'ects result in a final 0.5-parsec-long 
pillar for the isolated case whereas the final structure in the 
^^ciump/^^cioud=^ case is 1.5 parsec long. The driving process 
to form a pillar is how the matter distribution is able to first 
curve the shell and then to feed the base of the pillar by the hole 
mechanism identified in sect. 13. II (the holes are clearly visible in 
the final snapshots in Fig.fTOl). 
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Fig. 10: Density cuts of four simulations of an ionization front propagation on clumps (0.5 parsec in diameter). The UV flux at the 
top of the box is 10^ photons per second per square centimetre, the box is 4x2x2 parsecs (half of the cut is shown, only the central 
parsec of the simulation box). The last simulation (case(d)) is an isolated spherical clump of constant density 1000 cm^, the three 
others are clumps in a homogeneous medium of density 500 cm^. The clumps have densities of 1000 for case (a), 1750 for case (b) 
and 2500 cm^ for case (d). 
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Fig. 11: Monitoring of the size (top) and the mass (bottom) of 
pillars identified in Fig. [TOl The size is calculated by the diff'er- 
ence of the average position of the ionisation front and the its 
position at the centre of the y-z plane (head of the initial struc- 
ture). The mass is calculated in a box between the two vertical 
positions defined to calculated the size and with a width defined 
by the width of the initial structure and is normalized by the mass 
of the initial clump. 



The size evolution and the mass evolution of the pillars 
are comparable to the modulated-interface case. However the 
evolution is not linear due to the initial conditions. Indeed, 
there is a first phase in which the ionization front is curved and 
stretched verticaly around the clump when its propagation is 
slowed down by the overdensity. At this point (around 200 ky) 
the physical situation is comparable to the interface-modulated 
case, the ionization front is curved around a "hill". It also 
explains why we choose to do a longer simulation for the 
density-modulated case, it allows us to compare both cases at a 
"physically equivalent" state at the end. When the lateral shocks 
collide around the hill, the pillar captures a bit more than the 
mass of the initial clump and its mass increase slows down. This 
phase occurs around 300 ky and can be clearly identified in Fig. 
[TTl Then the mass increases due to the accumulation of matter 
at the base of the pillar. The process is less linear than in the 
interface-modulated case. 



A 



1 



1 
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Fig. 12: nHc/Mm/7/nHc/oMj=5 simulation at t=222 ky. Left: column 
density on a 0.4x0.4x0.4 pc^ box around the pillar structure in 
a face-on geometry. Right: mass-weighted histogram of the line 
of sight velocity in the same box (similar to optically-thin ob- 
servational line spectra). Each spectrum is made on a square of 
0.05x0.05 pc^ drawn on the column density map. The spectra are 
drawn between -4 and 4 km/s in 80 bins (horizontal axis) and the 
mass between 10"^ and 1 solar mass (vertical axis in log scale). 
The lateral shocks can be clearly identified in the double wings 
spectra leading to a very broad line width. 



4. Observational signature 

The two sections above present how pillars can be formed on a 
density or interface modulated region of the molecular cloud. 
These pillars are connected to the cloud, increasing in size 
and mass with time. However these variations can not be ob- 
served. A potentially observational signature to study is the 
str ucture of the line o f sight velocity. We use the same method 
as iGritschneder et al.l (l2010l) . We define a grid of squares of 
0.05x0.05 pc^ along the pillar and we bin the line of sight ve- 
locity in each square (see Fig. [121 [13] and [T4|. We use the simu- 
lation with the densest clump to study the line of sight velocity 
structure but the following results are generic. At t=222 ky (Fig. 
[12]), the lateral shocks can be clearly identified in the broad line 
spectrum. There are two components, a positive one coming to- 
wards us and a negative one going away. At this time the shell is 
curved around a hill which is comparable to the situation we got 
in the interface-modulated scenario after a short time. 

At t=444 ky (Fig. [13]), the lateral shocks have collided result- 
ing in a small line width for the histograms. The whole pillar is 
narrow and evacuates the pressure through radiative cooling and 
expansion to reach the equilibrium defined for the cold phase in 
figure [9] There is two dense part in the pillar: the head in which 
the matter of the clump has been accumulated and the base at 
the point where the lateral shocks have closed. It will resist the 
ionization because it is dense enough and matter is starting to 
accumulate at the base. 

At t= 998 ky (Fig. [14]), matter has accumulated at the base 
of the pillar, based on the process that was identified in Fig. 
[8] The line of sight velocity histograms have still a small line 
width compared to Fig[T2l This phenomenon is not specific 
to clumps simulations, the same analysis on the interface 
modulated simulation (aspect ratio of 0.9) in Fig. [15] and Fig. 
[T6] shows similar results but at an earlier stage. The interface 
modulated simulation presents double-wing velocity spectra 
around t=22.2 ky whereas the clumpy simulation presents these 
spectra around t=222 ky when the shock is curved around a hill. 
The double- wing structure of the velocity spectra is therefore a 
signature of the lateral shocks which are going to collide and a 
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Fig. 13: n}lciump/^^cioud=^ simulation at t=444 ky. Same plots 
as Fig.[T2l The lateral shocks collided and cannot been identified 
anymore. 

signature of the early stage of a forming pillar. 

The previous analysis can be applied to the diff'erent simula- 
tions presented in this paper without noticable change to our con- 
clusion. We also changed some physical conditions and param- 
eters. We added a constant external gravitational field, included 
self-gravity, doubled the resolution without global changes in the 
morphology of the final structures. We also changed the flux to 
an higher value (from 10^ to 5x10^ photons/s /cm^) and the struc- 
tures are similar but evolve faster as seen by lGritschneder et aP 
(120091) . The shell speed c shell is proportional to Fy^ (see Eq.[TT] 
and Eq. [121) therefore multiplying the flux by five increases the 
shell speed by 50% and the whole simulation evolves 50% faster. 
In all these situations, pillars form by the lateral collision of a 
curved shocked surface and the double- wing spectrum is always 
visible before the shocks collide. 

5. Conlusion and discussion 

We have presented a new scenario for the formation of structures 
at the edge of HIT regions and shown that: 

i A curved shock ahead of an ionization front can lead to a 
pillar if it is curved enough to collide lateraly on itself. 

ii This pocess is very efficient to form stable growing pillar, 
the narrower the initial structure, the more curved the front 
and the longer the pillar. 

iii Lateral gas flows can result in a density enhancement at the 
base of the pillars of 1-2 orders of magnitude compared to 
the collect and collapse scenario. 

iv When the shock is first formed flat, it can be curved 
by enough-dense clumps and leads to pillars. On isolated 
clumps (Radiation Driven Implosion) , the shock is naturally 
curved by the form of the clump but the resulting structure 
has a constant size and mass. 




Fig. 14: nllciump/^^cioud=^ simulation at t=998 ky. Same plots as 
Fig. [121 The line width is very small compared to the snapshots 
at t=222 ky. 

V The double- wing line spectra of the line-of- sight velocity is 
a signature of the lateral collision of the shock, and hence is a 
signature of the early stage of the formation of a pillar. It can 
be used as an observational signature for this new scenario. 

Various aspects of shocks orientation have been consid- 
ered by o ther authors. Oblique sho cks have previously been 
studied by Chevalier & Thevs (1975^ and cylindrical shocks by 
iKimura. Toshiya & Tosa (1990). These efl'ects have been fur- 
ther explored in the c ontext of 2D turbulent simulations (see 
lElmegreen et al.l [19951) . However these studies focus on the ef- 
fects of curvature on clumps enhancement. In this work, we have 
shown that shocks need to be enough curved to collide on them- 
selves in order to form a pillar like structure. 

At the edge of HII regions, the structure of the line-of- sight 
velocity can be investigated using radiotelescopes and suitable 
molecules tracing the dynamic, to detect the double- wing 
spectra on gas hills and thus nascent pillars. If a gas overdensity 
is detected at the top of the hill (e.g. using Herschel data), the 
shock curvature could be attributed to the presence of an initial 
clump, if not to the curvature of an initial interface. 

In this study, curved shocks have been generated using 
either density or interface modulations on a pc scale and in a 
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Fig. 15: w/h=0.9 simulation at t=22.2 ky (first column in Fig. 
[8]). The shock just formed on the modulation, the line width and 
the double- wing structure is comparable to Fig. [12] 
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Fig. 16: w/h=0.9 simulation at t= 133.2 ky (third column in Fig. 
(8]). The shock collided lateraly and the line width is small as in 
Fig. El 
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very simplified situation. However, we expect that the physical 
processes discussed in this paper can also be applied to more 
realistic situations where the ionization front interact with a tur- 
bulent interstellar medium. In this case the density fluctuations 
are much more complex, but the same mechanism should be 
at work to describe pillar formation. These turbulent situations 
will be studied in a forthcoming paper in which the interplay 
between shock curvature and turbulence will be studied to see 
its impact on star formation rates at the edge of HII regions. 



